# "Propaganda during economic crises: reference point adjustment in economic news" by Fatih Serkant Adiguzel


# This is the main script to produce all results except Figure 1.


rm(list = ls())

library(fixest)
library(dplyr)
library(magrittr)
library(here)
library(lubridate)
library(modelsummary)
library(ggplot2)
library(zoo)
library(purrr)
library(mgcv)
library(priceR)


load(here('main_data.rda'))

econ_news <- news %>% filter(econ_news == 1)

NONecon_news <- news %>% filter(econ_news == 0)


sabah <- news %>% filter(article_source == 'sabah') %>% filter(article_date < '2023-03-01')

sabah_summary <- sabah %>% group_by(article_date) %>% summarise(daily_articles = n())

# average
mean(sabah_summary$daily_articles)

trt <- news %>% filter(article_source == 'trt') %>% filter(article_date < '2023-03-01')

trt_summary <- trt %>% group_by(article_date) %>% summarise(daily_articles = n())

# average
mean(trt_summary$daily_articles)

sozcu <- news %>% filter(article_source == 'sozcu') %>% filter(article_date < '2023-03-01')
sozcu_summary <- sozcu %>% group_by(article_date) %>% summarise(daily_articles = n())

# average
mean(sozcu_summary$daily_articles)


# Create daily trends data

daily_trends_combined <- news %>% group_by(article_source, article_date) %>% summarise( world_news_count = sum(world_news, na.rm = TRUE),  world_news_count_wholetext = sum(world_news_wholetext, na.rm = TRUE), total_daily_articles = n(), foreign_econ_news_count = sum(foreign_econ_news, na.rm = TRUE), foreign_econ_news_count_wholetext = sum(foreign_econ_news_wholetext, na.rm = TRUE),
domestic_econ_news_count = sum(domestic_econ_news, na.rm = TRUE),
domestic_follows_foreignecon_count = sum(domestic_follows_foreignecon, na.rm = TRUE),
domestic_pos_follows_foreignecon_neg_count = sum(domestic_positive_follows_foreignecon_negative, na.rm = TRUE),
domestic_neut_follows_foreignecon_neg_count = sum(domestic_neutral_follows_foreignecon_negative, na.rm = TRUE))



daily_trends_combined$foreign_econ_news_share <- daily_trends_combined$foreign_econ_news_count / daily_trends_combined$total_daily_articles

daily_trends_combined$foreign_econ_news_share_wholetext <- daily_trends_combined$foreign_econ_news_count_wholetext / daily_trends_combined$total_daily_articles


daily_trends_combined$world_news_share <- daily_trends_combined$world_news_count / daily_trends_combined$total_daily_articles

daily_trends_combined$world_news_share_wholetext <- daily_trends_combined$world_news_count_wholetext / daily_trends_combined$total_daily_articles


daily_trends_combined$domestic_econ_news_share <- daily_trends_combined$domestic_econ_news_count / daily_trends_combined$total_daily_articles


daily_trends_combined$domestic_follows_foreignecon_share <- daily_trends_combined$domestic_follows_foreignecon_count / daily_trends_combined$total_daily_articles


daily_trends_combined$domestic_pos_neut_follows_foreign_neg_share <- (daily_trends_combined$domestic_pos_follows_foreignecon_neg_count + daily_trends_combined$domestic_neut_follows_foreignecon_neg_count) / daily_trends_combined$total_daily_articles



# Averages:
mean(daily_trends_combined$foreign_econ_news_share)
mean(daily_trends_combined$foreign_world_news_share)



daily_trends_combined$progov <- ifelse(daily_trends_combined$article_source == 'sozcu', 0 , 1)

daily_trends_combined_sabah_sozcu <- daily_trends_combined %>% filter(article_source %in% c('sozcu', 'sabah'))


daily_trends_combined_trt_sozcu <- daily_trends_combined %>% filter(article_source %in% c('sozcu', 'trt'))

#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

# Intl econ news coverage analysis (static) ----

econ_daily_trends_combined <- econ_news %>% group_by(article_source, article_date) %>% summarise(total_econ_daily_articles = n(), foreign_econ_news_count = sum(foreign_econ_news, na.rm = TRUE))

econ_daily_trends_combined$foreign_econ_news_share <- econ_daily_trends_combined$foreign_econ_news_count / econ_daily_trends_combined$total_econ_daily_articles


econ_daily_trends_combined$progov <- ifelse(econ_daily_trends_combined$article_source %in% c('sabah', 'trt'), 1, 0)


# Coverage models ----

coverage_all <- feols(foreign_econ_news_share~ progov | article_date, data=daily_trends_combined)

coverage_all_count <- feols(foreign_econ_news_count ~ progov | article_date, data=daily_trends_combined)


coverage_sabah <- feols(foreign_econ_news_share~ progov | article_date, data=daily_trends_combined_sabah_sozcu)

coverage_sabah_count <- feols(foreign_econ_news_count ~ progov | article_date, data=daily_trends_combined_sabah_sozcu)


coverage_trt <- feols(foreign_econ_news_share~ progov | article_date, data=daily_trends_combined_trt_sozcu)

coverage_trt_count <- feols(foreign_econ_news_count ~ progov | article_date, data=daily_trends_combined_trt_sozcu)


rows <- tribble(~term,  ~`1`,  ~`2`, ~`3`, ~`4`, ~`5`, ~`6`, 
                'Day FE', 'Yes', 'Yes', 'Yes', 'Yes',  'Yes', 'Yes',             'Average outcome', paste0(round(mean(daily_trends_combined$foreign_econ_news_share, na.rm = TRUE), digits = 2)), paste0(round(mean(daily_trends_combined_sabah_sozcu$foreign_econ_news_share, na.rm = TRUE), digits = 2)), paste0(round(mean(daily_trends_combined_trt_sozcu$foreign_econ_news_share, na.rm = TRUE), digits = 2)), paste0(round(mean(daily_trends_combined$foreign_econ_news_count, na.rm = TRUE), digits = 2)), paste0(round(mean(daily_trends_combined_sabah_sozcu$foreign_econ_news_count, na.rm = TRUE), digits = 2)), paste0(round(mean(daily_trends_combined_trt_sozcu$foreign_econ_news_count, na.rm = TRUE), digits = 2)))  


modelsummary(list(coverage_all, coverage_sabah, coverage_trt, coverage_all_count, coverage_sabah_count, coverage_trt_count), stars_note = TRUE, stars = c('*' = .1, '**' = .05, '***' = 0.01), metrics = 'R2', gof_omit = 'R2 Within|R2 Adj|Std.Errors|FE', add_rows = rows, coef_map = c('progov' = 'Pro-gov outlet'), output = here("table_1.txt")) 


#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

# Negative exposure models ----

# Calculate daily negative exposure in econ news first
international_econ_news <- news %>% filter(foreign_econ_news == 1)

daily_sentiment_trends <- international_econ_news %>% group_by(article_source, article_date) %>% summarise(total_econ_news_negativity = sum(negative_econ_news, na.rm = TRUE))

daily_trends_combined <- full_join(daily_trends_combined, daily_sentiment_trends, by = c('article_date', 'article_source'))

# If NA, then it means no negative exposure, so set it to 0.
daily_trends_combined$total_econ_news_negativity <- ifelse(is.na(daily_trends_combined$total_econ_news_negativity)==1, 0, daily_trends_combined$total_econ_news_negativity)


daily_trends_combined$total_negative_exposure <- daily_trends_combined$total_econ_news_negativity / daily_trends_combined$total_daily_articles

## define intl econ news wrt to the whole text this time:
international_econ_news_wholetext <- news %>% filter(foreign_econ_news_wholetext == 1)

daily_sentiment_trends_wholetext <- international_econ_news_wholetext %>% group_by(article_source, article_date) %>% summarise(total_econ_news_negativity_wholetext = sum(negative_econ_news, na.rm = TRUE))

daily_trends_combined <- full_join(daily_trends_combined, daily_sentiment_trends_wholetext, by = c('article_date', 'article_source'))


daily_trends_combined$total_econ_news_negativity_wholetext <- ifelse(is.na(daily_trends_combined$total_econ_news_negativity_wholetext)==1, 0, daily_trends_combined$total_econ_news_negativity_wholetext)


daily_trends_combined$total_negative_exposure_wholetext <- daily_trends_combined$total_econ_news_negativity_wholetext / daily_trends_combined$total_daily_articles



coverage_all_sent <- feols(total_negative_exposure ~ progov | article_date, data=daily_trends_combined)

coverage_all_count_sent <- feols(total_econ_news_negativity ~ progov | article_date, data=daily_trends_combined)

# re-run this due to newly added vars:
daily_trends_combined_sabah_sozcu <- daily_trends_combined %>% filter(article_source %in% c('sabah', 'sozcu'))

daily_trends_combined_trt_sozcu <- daily_trends_combined %>% filter(article_source %in% c('trt', 'sozcu'))


coverage_sabah_sent <- feols(total_negative_exposure ~ progov | article_date, data=daily_trends_combined_sabah_sozcu)

coverage_sabah_count_sent <- feols(total_econ_news_negativity ~ progov | article_date, data=daily_trends_combined_sabah_sozcu)


coverage_trt_sent <- feols(total_negative_exposure ~ progov | article_date, data=daily_trends_combined_trt_sozcu)

coverage_trt_count_sent <- feols(total_econ_news_negativity ~ progov | article_date, data=daily_trends_combined_trt_sozcu)



rows <- tribble(~term,  ~`1`,  ~`2`, ~`3`, ~`4`, ~`5`, ~`6`, 
     'Day FE', 'Yes', 'Yes', 'Yes', 'Yes',  'Yes', 'Yes',             'Average outcome', paste0(round(mean(daily_trends_combined$total_negative_exposure, na.rm = TRUE), digits = 2)), paste0(round(mean(daily_trends_combined_sabah_sozcu$total_negative_exposure, na.rm = TRUE), digits = 2)), paste0(round(mean(daily_trends_combined_trt_sozcu$total_negative_exposure, na.rm = TRUE), digits = 2)), paste0(round(mean(daily_trends_combined$total_econ_news_negativity, na.rm = TRUE), digits = 2)), paste0(round(mean(daily_trends_combined_sabah_sozcu$total_econ_news_negativity, na.rm = TRUE), digits = 2)), paste0(round(mean(daily_trends_combined_trt_sozcu$total_econ_news_negativity, na.rm = TRUE), digits = 2)))  

modelsummary(list(coverage_all_sent, coverage_sabah_sent, coverage_trt_sent, coverage_all_count_sent, coverage_sabah_count_sent, coverage_trt_count_sent), stars_note = TRUE, stars = c('*' = .1, '**' = .05, '***' = 0.01), metrics = 'R2', gof_omit = 'R2 Within|R2 Adj|Std.Errors|FE', add_rows = rows, coef_map = c('progov' = 'Pro-gov outlet'), output = here('table_2.txt')) 


# Descriptive Figures ----


econ_coverage_plot <- ggplot(daily_trends_combined, aes(x = article_date, y = foreign_econ_news_share, color = article_source)) + geom_point(size = 0.4, alpha = 0.7) + geom_smooth(method = 'loess', span = 0.1, se = FALSE) + geom_vline(xintercept =c(as.Date('2021-09-24'), as.Date('2021-10-22'), as.Date('2021-09-24'), as.Date('2021-10-22'), as.Date('2021-11-19'), as.Date('2021-12-17'), as.Date('2022-08-19'), as.Date('2022-09-23'), as.Date('2022-10-21'), as.Date('2022-11-25'), as.Date('2023-02-23'))) + theme_bw() + theme(legend.position = 'bottom', axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + labs(y = 'Daily foreign economy coverage', x = '', color = 'newspaper') + scale_x_date(date_breaks = 'month') 


ggsave(econ_coverage_plot, filename = here('figure_2.pdf'), width = 7, height = 6)



econ_sentiment_plot <- ggplot(daily_trends_combined, aes(x = article_date, y = total_negative_exposure, color = factor(article_source))) + geom_point(size = 0.4, alpha = 0.7) + geom_smooth(method = 'loess', span = 0.1, se = FALSE) + geom_vline(xintercept =c(as.Date('2021-09-24'), as.Date('2021-10-22'), as.Date('2021-09-24'), as.Date('2021-10-22'), as.Date('2021-11-19'), as.Date('2021-12-17'), as.Date('2022-08-19'), as.Date('2022-09-23'), as.Date('2022-10-21'), as.Date('2022-11-25'), as.Date('2023-02-23')))+ theme_bw() + theme(legend.position = 'bottom', axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + labs(y = 'Daily negative sentiment exposure in foreign economy news', x = '', color = 'newspaper') + scale_x_date(date_breaks = 'month') 

ggsave(econ_sentiment_plot, filename = here('figure_3.pdf'), width = 7, height = 6)


# Descriptive Figures (Appendix) ----

#### Currency plot ----

d <- read.csv(here('currency_data.csv'))

d$date <- as.Date(d$date)

currency_plot <- ggplot(data = d, aes(x = date, y = usdtry))+geom_point(size = 0.05)+geom_line()+ scale_x_date(breaks = '3 months')+ theme_classic()+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+ylab('USD/TRY')+xlab('')+geom_vline(xintercept = as.Date('2021-08-04'), linetype = 'dashed', color ='red') + geom_vline(xintercept = as.Date('2021-09-23'), linetype = 'dashed', color ='black')

ggsave(currency_plot, filename = here('figure_A1.pdf'), width = 7, height = 6)

#### Coverage plot ----

coverage_trend_plot <- ggplot(daily_trends_combined, aes(x = article_date, y = total_daily_articles, color = article_source)) + geom_point(size = 0.4, alpha = 0.7) + geom_smooth(method = 'loess', span = 0.3, se = FALSE) + geom_vline(xintercept =c(as.Date('2021-09-24'), as.Date('2021-10-22'), as.Date('2021-09-24'), as.Date('2021-10-22'), as.Date('2021-11-19'), as.Date('2021-12-17'), as.Date('2022-08-19'), as.Date('2022-09-23'), as.Date('2022-10-21'), as.Date('2022-11-25'), as.Date('2023-02-23'))) + theme_bw() + theme(legend.position = 'bottom', axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + labs(y = 'Total number of articles (daily)', x = '', color = 'newspaper') + scale_x_date(date_breaks = 'month') 

ggsave(coverage_trend_plot, file = here('figure_A2.pdf'), width = 7, height = 6)


#### Domestic econ coverage plot ----

domestic_econ_coverage_plot <- ggplot(daily_trends_combined, aes(x = article_date, y = domestic_econ_news_share, color = article_source)) + geom_point(size = 0.4, alpha = 0.7) + geom_smooth(method = 'loess', span = 0.1, se = FALSE) + geom_vline(xintercept =c(as.Date('2021-09-24'), as.Date('2021-10-22'), as.Date('2021-09-24'), as.Date('2021-10-22'), as.Date('2021-11-19'), as.Date('2021-12-17'), as.Date('2022-08-19'), as.Date('2022-09-23'), as.Date('2022-10-21'), as.Date('2022-11-25'), as.Date('2023-02-23'))) + theme_bw() + theme(legend.position = 'bottom', axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + labs(y = 'Daily domestic economy coverage', x = '', color = 'newspaper') + scale_x_date(date_breaks = 'month') 

ggsave(domestic_econ_coverage_plot, filename = here('figure_A6.pdf'), width = 7, height = 6)

#### World news coverage plot ----

# note that there was an issue with TRT data on December 31, 2022  and there is only 1 article scraped for that day, which is classified as foreign news. I got rid of that one observation below since it changes the scale dramatically

world_news_coverage_plot <- ggplot(subset(daily_trends_combined,world_news_share!=1), aes(x = article_date, y = world_news_share, color = article_source)) + geom_point(size = 0.4, alpha = 0.7) + geom_smooth(method = 'loess', span = 0.1, se = FALSE) + geom_vline(xintercept =c(as.Date('2021-09-24'), as.Date('2021-10-22'), as.Date('2021-09-24'), as.Date('2021-10-22'), as.Date('2021-11-19'), as.Date('2021-12-17'), as.Date('2022-08-19'), as.Date('2022-09-23'), as.Date('2022-10-21'), as.Date('2022-11-25'), as.Date('2023-02-23'))) + theme_bw() + theme(legend.position = 'bottom', axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + labs(y = 'Daily world news coverage (excluding foreign economy news)', x = '', color = 'newspaper') + scale_x_date(date_breaks = 'month') 

ggsave(world_news_coverage_plot, filename = here('figure_B1.pdf'), width = 7, height = 6)


### Topic Distribution over time ----

#### Sabah topics ----


daily_article_trend_sabah <- sabah %>% group_by(article_date) %>% summarise(total_daily_econ = sum(econ_news, na.rm = TRUE), total_daily_culture = sum(culture_news, na.rm = TRUE), total_daily_health = sum(health_news, na.rm = TRUE), total_daily_politics = sum(politics_news, na.rm = TRUE), total_daily_sports = sum(sports_news, na.rm = TRUE), total_daily_technology = sum(tech_news, na.rm = TRUE), total_daily_world = sum(world_news, na.rm = TRUE), total_daily_education = sum(education_news, na.rm = TRUE), total_daily_local = sum(local_news, na.rm = TRUE), total_daily_magazine = sum(magazine_news, na.rm = TRUE), total_daily_other = sum(other_news, na.rm = TRUE), total_daily_weather_disaster = sum(weather_disaster_news, na.rm = TRUE), total_daily_articles = n()) %>% ungroup()

daily_article_trend_sabah$econtopic_daily_share <- daily_article_trend_sabah$total_daily_econ / daily_article_trend_sabah$total_daily_articles

daily_article_trend_sabah$politicstopic_daily_share <- daily_article_trend_sabah$total_daily_politics / daily_article_trend_sabah$total_daily_articles


daily_article_trend_sabah$culturetopic_daily_share <- daily_article_trend_sabah$total_daily_culture / daily_article_trend_sabah$total_daily_articles


daily_article_trend_sabah$healthtopic_daily_share <- daily_article_trend_sabah$total_daily_health / daily_article_trend_sabah$total_daily_articles

daily_article_trend_sabah$sporttopic_daily_share <- daily_article_trend_sabah$total_daily_sports / daily_article_trend_sabah$total_daily_articles

daily_article_trend_sabah$technologytopic_daily_share <- daily_article_trend_sabah$total_daily_technology / daily_article_trend_sabah$total_daily_articles

daily_article_trend_sabah$worldtopic_daily_share <- daily_article_trend_sabah$total_daily_world / daily_article_trend_sabah$total_daily_articles


daily_article_trend_sabah$educationtopic_daily_share <- daily_article_trend_sabah$total_daily_education / daily_article_trend_sabah$total_daily_articles

daily_article_trend_sabah$localtopic_daily_share <- daily_article_trend_sabah$total_daily_local / daily_article_trend_sabah$total_daily_articles

daily_article_trend_sabah$magazinetopic_daily_share <- daily_article_trend_sabah$total_daily_magazine / daily_article_trend_sabah$total_daily_articles

daily_article_trend_sabah$othertopic_daily_share <- daily_article_trend_sabah$total_daily_other / daily_article_trend_sabah$total_daily_articles

daily_article_trend_sabah$weatherdisastertopic_daily_share <- daily_article_trend_sabah$total_daily_weather_disaster / daily_article_trend_sabah$total_daily_articles


daily_article_trend_sabah_long <- daily_article_trend_sabah %>% select(article_date, econtopic_daily_share, culturetopic_daily_share, educationtopic_daily_share, localtopic_daily_share, magazinetopic_daily_share, othertopic_daily_share,  healthtopic_daily_share, politicstopic_daily_share, worldtopic_daily_share, sporttopic_daily_share, technologytopic_daily_share, weatherdisastertopic_daily_share)


daily_article_trend_sabah_long <- daily_article_trend_sabah_long %>% tidyr::pivot_longer(econtopic_daily_share:weatherdisastertopic_daily_share) %>% unique() 

daily_article_trend_sabah_long$topics <- ifelse(daily_article_trend_sabah_long$name=='culturetopic_daily_share', 'Culture', ifelse(daily_article_trend_sabah_long$name=='econtopic_daily_share', 'Economy', ifelse(daily_article_trend_sabah_long$name=='educationtopic_daily_share', 'Education', ifelse(daily_article_trend_sabah_long$name=='healthtopic_daily_share', 'Health', ifelse(daily_article_trend_sabah_long$name=='localtopic_daily_share', 'Local', ifelse(daily_article_trend_sabah_long$name == 'magazinetopic_daily_share', 'Magazine', ifelse(daily_article_trend_sabah_long$name == 'othertopic_daily_share', 'Other', ifelse(daily_article_trend_sabah_long$name == 'politicstopic_daily_share', 'Politics', ifelse(daily_article_trend_sabah_long$name == 'sporttopic_daily_share', 'Sports', ifelse(daily_article_trend_sabah_long$name == 'technologytopic_daily_share', 'Technology', ifelse(daily_article_trend_sabah_long$name == 'weatherdisastertopic_daily_share', 'Weather/disaster', ifelse(daily_article_trend_sabah_long$name == 'worldtopic_daily_share', 'World', NA))))))))))))


all_topics_panel_daily_sabah <- ggplot(data = daily_article_trend_sabah_long, aes(x = article_date, y = value, color = as.factor(topics))) + geom_point(size = 0.03, alpha = 0.2) + geom_smooth(method = "loess", span = 0.3) + scale_x_date(breaks = '3 months') + facet_wrap(~ topics, scales = 'free') + theme_bw() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5), legend.position = "none") + labs(y = 'Topic share (daily)', x = '', color = 'Topics')

ggsave(all_topics_panel_daily_sabah, file = here('figure_A3.pdf'), width = 8, height = 8)


#### Sozcu topics ----

daily_article_trend_sozcu <- sozcu %>% group_by(article_date) %>% summarise(total_daily_econ = sum(econ_news, na.rm = TRUE), total_daily_culture = sum(culture_news, na.rm = TRUE), total_daily_health = sum(health_news, na.rm = TRUE), total_daily_politics = sum(politics_news, na.rm = TRUE), total_daily_sports = sum(sports_news, na.rm = TRUE), total_daily_technology = sum(tech_news, na.rm = TRUE), total_daily_world = sum(world_news, na.rm = TRUE), total_daily_education = sum(education_news, na.rm = TRUE), total_daily_local = sum(local_news, na.rm = TRUE), total_daily_magazine = sum(magazine_news, na.rm = TRUE), total_daily_other = sum(other_news, na.rm = TRUE), total_daily_weather_disaster = sum(weather_disaster_news, na.rm = TRUE), total_daily_articles = n()) %>% ungroup()

daily_article_trend_sozcu$econtopic_daily_share <- daily_article_trend_sozcu$total_daily_econ / daily_article_trend_sozcu$total_daily_articles

daily_article_trend_sozcu$politicstopic_daily_share <- daily_article_trend_sozcu$total_daily_politics / daily_article_trend_sozcu$total_daily_articles


daily_article_trend_sozcu$culturetopic_daily_share <- daily_article_trend_sozcu$total_daily_culture / daily_article_trend_sozcu$total_daily_articles


daily_article_trend_sozcu$healthtopic_daily_share <- daily_article_trend_sozcu$total_daily_health / daily_article_trend_sozcu$total_daily_articles

daily_article_trend_sozcu$sporttopic_daily_share <- daily_article_trend_sozcu$total_daily_sports / daily_article_trend_sozcu$total_daily_articles

daily_article_trend_sozcu$technologytopic_daily_share <- daily_article_trend_sozcu$total_daily_technology / daily_article_trend_sozcu$total_daily_articles

daily_article_trend_sozcu$worldtopic_daily_share <- daily_article_trend_sozcu$total_daily_world / daily_article_trend_sozcu$total_daily_articles


daily_article_trend_sozcu$educationtopic_daily_share <- daily_article_trend_sozcu$total_daily_education / daily_article_trend_sozcu$total_daily_articles

daily_article_trend_sozcu$localtopic_daily_share <- daily_article_trend_sozcu$total_daily_local / daily_article_trend_sozcu$total_daily_articles

daily_article_trend_sozcu$magazinetopic_daily_share <- daily_article_trend_sozcu$total_daily_magazine / daily_article_trend_sozcu$total_daily_articles

daily_article_trend_sozcu$othertopic_daily_share <- daily_article_trend_sozcu$total_daily_other / daily_article_trend_sozcu$total_daily_articles

daily_article_trend_sozcu$weatherdisastertopic_daily_share <- daily_article_trend_sozcu$total_daily_weather_disaster / daily_article_trend_sozcu$total_daily_articles


daily_article_trend_sozcu_long <- daily_article_trend_sozcu %>% select(article_date, econtopic_daily_share, culturetopic_daily_share, educationtopic_daily_share, localtopic_daily_share, magazinetopic_daily_share, othertopic_daily_share,  healthtopic_daily_share, politicstopic_daily_share, worldtopic_daily_share, sporttopic_daily_share, technologytopic_daily_share, weatherdisastertopic_daily_share)


daily_article_trend_sozcu_long <- daily_article_trend_sozcu_long %>% tidyr::pivot_longer(econtopic_daily_share:weatherdisastertopic_daily_share) %>% unique() 

daily_article_trend_sozcu_long$topics <- ifelse(daily_article_trend_sozcu_long$name=='culturetopic_daily_share', 'Culture', ifelse(daily_article_trend_sozcu_long$name=='econtopic_daily_share', 'Economy', ifelse(daily_article_trend_sozcu_long$name=='educationtopic_daily_share', 'Education', ifelse(daily_article_trend_sozcu_long$name=='healthtopic_daily_share', 'Health', ifelse(daily_article_trend_sozcu_long$name=='localtopic_daily_share', 'Local', ifelse(daily_article_trend_sozcu_long$name == 'magazinetopic_daily_share', 'Magazine', ifelse(daily_article_trend_sozcu_long$name == 'othertopic_daily_share', 'Other', ifelse(daily_article_trend_sozcu_long$name == 'politicstopic_daily_share', 'Politics', ifelse(daily_article_trend_sozcu_long$name == 'sporttopic_daily_share', 'Sports', ifelse(daily_article_trend_sozcu_long$name == 'technologytopic_daily_share', 'Technology', ifelse(daily_article_trend_sozcu_long$name == 'weatherdisastertopic_daily_share', 'Weather/disaster', ifelse(daily_article_trend_sozcu_long$name == 'worldtopic_daily_share', 'World', NA))))))))))))


all_topics_panel_daily_sozcu <- ggplot(data = daily_article_trend_sozcu_long, aes(x = article_date, y = value, color = as.factor(topics))) + geom_point(size = 0.03, alpha = 0.2) + geom_smooth(method = "loess", span = 0.3) + scale_x_date(breaks = '3 months') + facet_wrap(~ topics, scales = 'free') + theme_bw() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5), legend.position = "none") + labs(y = 'Topic share (daily)', x = '', color = 'Topics')

ggsave(all_topics_panel_daily_sozcu, file = here('figure_A5.pdf'), width = 8, height = 8)



#### TRT topics ----
daily_article_trend_trt <- trt %>% group_by(article_date) %>% summarise(total_daily_econ = sum(econ_news, na.rm = TRUE), total_daily_culture = sum(culture_news, na.rm = TRUE), total_daily_health = sum(health_news, na.rm = TRUE), total_daily_politics = sum(politics_news, na.rm = TRUE), total_daily_sports = sum(sports_news, na.rm = TRUE), total_daily_technology = sum(tech_news, na.rm = TRUE), total_daily_world = sum(world_news, na.rm = TRUE), total_daily_education = sum(education_news, na.rm = TRUE), total_daily_local = sum(local_news, na.rm = TRUE), total_daily_magazine = sum(magazine_news, na.rm = TRUE), total_daily_other = sum(other_news, na.rm = TRUE), total_daily_weather_disaster = sum(weather_disaster_news, na.rm = TRUE), total_daily_articles = n()) %>% ungroup()

daily_article_trend_trt$econtopic_daily_share <- daily_article_trend_trt$total_daily_econ / daily_article_trend_trt$total_daily_articles

daily_article_trend_trt$politicstopic_daily_share <- daily_article_trend_trt$total_daily_politics / daily_article_trend_trt$total_daily_articles


daily_article_trend_trt$culturetopic_daily_share <- daily_article_trend_trt$total_daily_culture / daily_article_trend_trt$total_daily_articles


daily_article_trend_trt$healthtopic_daily_share <- daily_article_trend_trt$total_daily_health / daily_article_trend_trt$total_daily_articles

daily_article_trend_trt$sporttopic_daily_share <- daily_article_trend_trt$total_daily_sports / daily_article_trend_trt$total_daily_articles

daily_article_trend_trt$technologytopic_daily_share <- daily_article_trend_trt$total_daily_technology / daily_article_trend_trt$total_daily_articles

daily_article_trend_trt$worldtopic_daily_share <- daily_article_trend_trt$total_daily_world / daily_article_trend_trt$total_daily_articles


daily_article_trend_trt$educationtopic_daily_share <- daily_article_trend_trt$total_daily_education / daily_article_trend_trt$total_daily_articles

daily_article_trend_trt$localtopic_daily_share <- daily_article_trend_trt$total_daily_local / daily_article_trend_trt$total_daily_articles

daily_article_trend_trt$magazinetopic_daily_share <- daily_article_trend_trt$total_daily_magazine / daily_article_trend_trt$total_daily_articles

daily_article_trend_trt$othertopic_daily_share <- daily_article_trend_trt$total_daily_other / daily_article_trend_trt$total_daily_articles

daily_article_trend_trt$weatherdisastertopic_daily_share <- daily_article_trend_trt$total_daily_weather_disaster / daily_article_trend_trt$total_daily_articles


daily_article_trend_trt_long <- daily_article_trend_trt %>% select(article_date, econtopic_daily_share, culturetopic_daily_share, educationtopic_daily_share, localtopic_daily_share, magazinetopic_daily_share, othertopic_daily_share,  healthtopic_daily_share, politicstopic_daily_share, worldtopic_daily_share, sporttopic_daily_share, technologytopic_daily_share, weatherdisastertopic_daily_share)


daily_article_trend_trt_long <- daily_article_trend_trt_long %>% tidyr::pivot_longer(econtopic_daily_share:weatherdisastertopic_daily_share) %>% unique() 

daily_article_trend_trt_long$topics <- ifelse(daily_article_trend_trt_long$name=='culturetopic_daily_share', 'Culture', ifelse(daily_article_trend_trt_long$name=='econtopic_daily_share', 'Economy', ifelse(daily_article_trend_trt_long$name=='educationtopic_daily_share', 'Education', ifelse(daily_article_trend_trt_long$name=='healthtopic_daily_share', 'Health', ifelse(daily_article_trend_trt_long$name=='localtopic_daily_share', 'Local', ifelse(daily_article_trend_trt_long$name == 'magazinetopic_daily_share', 'Magazine', ifelse(daily_article_trend_trt_long$name == 'othertopic_daily_share', 'Other', ifelse(daily_article_trend_trt_long$name == 'politicstopic_daily_share', 'Politics', ifelse(daily_article_trend_trt_long$name == 'sporttopic_daily_share', 'Sports', ifelse(daily_article_trend_trt_long$name == 'technologytopic_daily_share', 'Technology', ifelse(daily_article_trend_trt_long$name == 'weatherdisastertopic_daily_share', 'Weather/disaster', ifelse(daily_article_trend_trt_long$name == 'worldtopic_daily_share', 'World', NA))))))))))))


all_topics_panel_daily_trt <- ggplot(data = daily_article_trend_trt_long, aes(x = article_date, y = value, color = as.factor(topics))) + geom_point(size = 0.03, alpha = 0.2) + geom_smooth(method = "loess", span = 0.3) + scale_x_date(breaks = '3 months') + facet_wrap(~ topics, scales = 'free') + theme_bw() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5), legend.position = "none") + labs(y = 'Topic share (daily)', x = '', color = 'Topics')

ggsave(all_topics_panel_daily_trt, file = here('Figure_A4.pdf'), width = 8, height = 8)



#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

# Foreign economy coverage analysis ----

daily_trends_combined$year <- year(daily_trends_combined$article_date)


daily_trends_combined$month <- month(daily_trends_combined$article_date)

daily_trends_combined$month_year <- paste0(daily_trends_combined$year,'-', daily_trends_combined$month)

daily_trends_combined$month_yearDate <- as.yearmon(daily_trends_combined$month_year)

daily_trends_combined$month_yearDate <- as.Date(daily_trends_combined$month_yearDate)



#### Main: (baseline = 2021 December) ----


# Note that all dynamic models in this script use day fixed effects. Therefore, monthly dummies will be dropped since models only exploit within-day variation. Warning is expected here.

coverage_foreignecon_dynamic <- feols(foreign_econ_news_share~ progov*i(month_yearDate, as.Date('2021-12-01')) | article_date, data=daily_trends_combined, panel.id = c('article_source', 'article_date'), vcov = 'DK')

# Extract coefs for the plot:
coverage_foreignecon_dynamic_result <- summary(coverage_foreignecon_dynamic)

coverage_foreignecon_dynamic_result_data <- data.frame(coefs = coverage_foreignecon_dynamic_result$coefficients, se = coverage_foreignecon_dynamic_result$se) %>% tibble::rownames_to_column() %>% filter(grepl('progov:month_yearDate', rowname))


coverage_foreignecon_dynamic_result_data$lo95 <- coverage_foreignecon_dynamic_result_data$coefs - 1.96*coverage_foreignecon_dynamic_result_data$se 

coverage_foreignecon_dynamic_result_data$hi95 <- coverage_foreignecon_dynamic_result_data$coefs + 1.96*coverage_foreignecon_dynamic_result_data$se 

coverage_foreignecon_dynamic_result_data$variable_name <- stringr::str_split_fixed(coverage_foreignecon_dynamic_result_data$rowname, '::', 3)[, 2]

coverage_foreignecon_dynamic_result_data$variable_name <- as.Date(coverage_foreignecon_dynamic_result_data$variable_name)

coverage_foreignecon_dynamic_result_data <- coverage_foreignecon_dynamic_result_data %>% add_row(rowname = 'baseline', coefs = 0, se = 0, lo95 = 0, hi95 = 0, variable_name = as.Date('2021-12-01'))


foreign_econ_dynamic_plot <- ggplot(coverage_foreignecon_dynamic_result_data, aes(x=variable_name, y = coefs, ymin = lo95, ymax = hi95))+geom_point()+geom_linerange()+theme_bw()+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+ylab('Estimated effect of publishing an article on foreign economy')+ xlab('')+geom_hline(yintercept = 0, color = 'red', linetype = 'dashed')+ geom_rect(xmin = as.Date('2021-09-24'), xmax = as.Date('2021-12-17'), ymin=-Inf, ymax=Inf, alpha=0.005, fill="firebrick") + geom_rect(xmin = as.Date('2022-08-19'), xmax = as.Date('2022-11-25'), ymin=-Inf, ymax=Inf, alpha=0.005, fill="firebrick")+ scale_x_date(date_breaks = "1 month", date_labels = "%Y-%m")

ggsave(foreign_econ_dynamic_plot, file = here('figure_4.pdf'), width = 7, height = 6)

## substantive interpretation
coverage_foreignecon_dynamic_result_data[12:23, ] %>% pull(coefs) %>% mean()


### Sabah vs Sozcu (baseline = 2021 December) ----

daily_trends_sabahsozcu <- daily_trends_combined %>% filter(article_source!='trt')

coverage_dynamic_sabahsozcu <- feols(foreign_econ_news_share~ progov*i(month_yearDate, ref = as.Date('2021-12-01')) | article_date, data=daily_trends_sabahsozcu, panel.id = c('article_source', 'article_date'), vcov = 'DK')


coverage_dynamic_sabahsozcu_data <- data.frame(coefs = coverage_dynamic_sabahsozcu$coefficients, se = coverage_dynamic_sabahsozcu$se) %>% tibble::rownames_to_column() %>% filter(grepl('progov:month_yearDate', rowname))


coverage_dynamic_sabahsozcu_data$lo95 <- coverage_dynamic_sabahsozcu_data$coefs - 1.96*coverage_dynamic_sabahsozcu_data$se 

coverage_dynamic_sabahsozcu_data$hi95 <- coverage_dynamic_sabahsozcu_data$coefs + 1.96*coverage_dynamic_sabahsozcu_data$se 

coverage_dynamic_sabahsozcu_data$variable_name <- stringr::str_split_fixed(coverage_dynamic_sabahsozcu_data$rowname, '::', 3)[, 2]
coverage_dynamic_sabahsozcu_data$variable_name <- as.Date(coverage_dynamic_sabahsozcu_data$variable_name)

coverage_dynamic_sabahsozcu_data <- coverage_dynamic_sabahsozcu_data %>% add_row(rowname = 'baseline', coefs = 0, se = 0, lo95 = 0, hi95 = 0, variable_name = as.Date('2021-12-01'))


coverage_dynamic_sabahsozcu_plot <- ggplot(coverage_dynamic_sabahsozcu_data, aes(x=variable_name, y = coefs, ymin = lo95, ymax = hi95))+geom_point()+geom_linerange()+theme_bw()+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+ylab('Estimated effect of publishing an article on foreign economy')+ xlab('')+geom_hline(yintercept = 0, color = 'red', linetype = 'dashed')+ geom_rect(xmin = as.Date('2021-09-24'), xmax = as.Date('2021-12-17'), ymin=-Inf, ymax=Inf, alpha=0.005, fill="firebrick") + geom_rect(xmin = as.Date('2022-08-19'), xmax = as.Date('2022-11-25'), ymin=-Inf, ymax=Inf, alpha=0.005, fill="firebrick")+ scale_x_date(date_breaks = "1 month", date_labels = "%Y-%m")

ggsave(coverage_dynamic_sabahsozcu_plot, file = here('figure_C1.pdf'), width = 7, height = 6)


### TRT vs Sozcu (baseline = 2021 December) ----

daily_trends_trtsozcu <- daily_trends_combined %>% filter(article_source!='sabah')

coverage_dynamic_trtsozcu <- feols(foreign_econ_news_share~ progov*i(month_yearDate, ref = as.Date('2021-12-01')) | article_date, data=daily_trends_trtsozcu, panel.id = c('article_source', 'article_date'), vcov = 'DK')


coverage_dynamic_trtsozcu_data <- data.frame(coefs = coverage_dynamic_trtsozcu$coefficients, se = coverage_dynamic_trtsozcu$se) %>% tibble::rownames_to_column() %>% filter(grepl('progov:month_yearDate', rowname))


coverage_dynamic_trtsozcu_data$lo95 <- coverage_dynamic_trtsozcu_data$coefs - 1.96*coverage_dynamic_trtsozcu_data$se 

coverage_dynamic_trtsozcu_data$hi95 <- coverage_dynamic_trtsozcu_data$coefs + 1.96*coverage_dynamic_trtsozcu_data$se 

coverage_dynamic_trtsozcu_data$variable_name <- stringr::str_split_fixed(coverage_dynamic_trtsozcu_data$rowname, '::', 3)[, 2]
coverage_dynamic_trtsozcu_data$variable_name <- as.Date(coverage_dynamic_trtsozcu_data$variable_name)

coverage_dynamic_trtsozcu_data <- coverage_dynamic_trtsozcu_data %>% add_row(rowname = 'baseline', coefs = 0, se = 0, lo95 = 0, hi95 = 0, variable_name = as.Date('2021-12-01'))


coverage_dynamic_trtsozcu_plot <- ggplot(coverage_dynamic_trtsozcu_data, aes(x=variable_name, y = coefs, ymin = lo95, ymax = hi95))+geom_point()+geom_linerange()+theme_bw()+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+ylab('Estimated effect of publishing an article on foreign economy')+ xlab('')+geom_hline(yintercept = 0, color = 'red', linetype = 'dashed')+ geom_rect(xmin = as.Date('2021-09-24'), xmax = as.Date('2021-12-17'), ymin=-Inf, ymax=Inf, alpha=0.005, fill="firebrick") + geom_rect(xmin = as.Date('2022-08-19'), xmax = as.Date('2022-11-25'), ymin=-Inf, ymax=Inf, alpha=0.005, fill="firebrick")+ scale_x_date(date_breaks = "1 month", date_labels = "%Y-%m")

ggsave(coverage_dynamic_trtsozcu_plot, file = here('figure_C2.pdf'), width = 7, height = 6)


#### Robustness: (baseline = 2021 November) ----

coverage_foreignecon_dynamic <- feols(foreign_econ_news_share~ progov*i(month_yearDate, as.Date('2021-11-01')) | article_date, data=daily_trends_combined, panel.id = c('article_source', 'article_date'), vcov = 'DK')

coverage_foreignecon_dynamic_result <- summary(coverage_foreignecon_dynamic)

coverage_foreignecon_dynamic_result_data <- data.frame(coefs = coverage_foreignecon_dynamic_result$coefficients, se = coverage_foreignecon_dynamic_result$se) %>% tibble::rownames_to_column() %>% filter(grepl('progov:month_yearDate', rowname))


coverage_foreignecon_dynamic_result_data$lo95 <- coverage_foreignecon_dynamic_result_data$coefs - 1.96*coverage_foreignecon_dynamic_result_data$se 

coverage_foreignecon_dynamic_result_data$hi95 <- coverage_foreignecon_dynamic_result_data$coefs + 1.96*coverage_foreignecon_dynamic_result_data$se 

coverage_foreignecon_dynamic_result_data$variable_name <- stringr::str_split_fixed(coverage_foreignecon_dynamic_result_data$rowname, '::', 3)[, 2]

coverage_foreignecon_dynamic_result_data$variable_name <- as.Date(coverage_foreignecon_dynamic_result_data$variable_name)

coverage_foreignecon_dynamic_result_data <- coverage_foreignecon_dynamic_result_data %>% add_row(rowname = 'baseline', coefs = 0, se = 0, lo95 = 0, hi95 = 0, variable_name = as.Date('2021-11-01'))


foreign_econ_dynamic_plot <- ggplot(coverage_foreignecon_dynamic_result_data, aes(x=variable_name, y = coefs, ymin = lo95, ymax = hi95))+geom_point()+geom_linerange()+theme_bw()+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+ylab('Estimated effect of publishing an article on foreign economy')+ xlab('')+geom_hline(yintercept = 0, color = 'red', linetype = 'dashed')+ geom_rect(xmin = as.Date('2021-09-24'), xmax = as.Date('2021-12-17'), ymin=-Inf, ymax=Inf, alpha=0.005, fill="firebrick") + geom_rect(xmin = as.Date('2022-08-19'), xmax = as.Date('2022-11-25'), ymin=-Inf, ymax=Inf, alpha=0.005, fill="firebrick")+ scale_x_date(date_breaks = "1 month", date_labels = "%Y-%m")

ggsave(foreign_econ_dynamic_plot, file = here('figure_C3.pdf'), width = 7, height = 6)


#### Robustness: (baseline = 2021 October) ----

coverage_foreignecon_dynamic <- feols(foreign_econ_news_share~ progov*i(month_yearDate, as.Date('2021-10-01')) | article_date, data=daily_trends_combined, panel.id = c('article_source', 'article_date'), vcov = 'DK')

coverage_foreignecon_dynamic_result <- summary(coverage_foreignecon_dynamic)

coverage_foreignecon_dynamic_result_data <- data.frame(coefs = coverage_foreignecon_dynamic_result$coefficients, se = coverage_foreignecon_dynamic_result$se) %>% tibble::rownames_to_column() %>% filter(grepl('progov:month_yearDate', rowname))


coverage_foreignecon_dynamic_result_data$lo95 <- coverage_foreignecon_dynamic_result_data$coefs - 1.96*coverage_foreignecon_dynamic_result_data$se 

coverage_foreignecon_dynamic_result_data$hi95 <- coverage_foreignecon_dynamic_result_data$coefs + 1.96*coverage_foreignecon_dynamic_result_data$se 

coverage_foreignecon_dynamic_result_data$variable_name <- stringr::str_split_fixed(coverage_foreignecon_dynamic_result_data$rowname, '::', 3)[, 2]

coverage_foreignecon_dynamic_result_data$variable_name <- as.Date(coverage_foreignecon_dynamic_result_data$variable_name)

coverage_foreignecon_dynamic_result_data <- coverage_foreignecon_dynamic_result_data %>% add_row(rowname = 'baseline', coefs = 0, se = 0, lo95 = 0, hi95 = 0, variable_name = as.Date('2021-10-01'))


foreign_econ_dynamic_plot <- ggplot(coverage_foreignecon_dynamic_result_data, aes(x=variable_name, y = coefs, ymin = lo95, ymax = hi95))+geom_point()+geom_linerange()+theme_bw()+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+ylab('Estimated effect of publishing an article on foreign economy')+ xlab('')+geom_hline(yintercept = 0, color = 'red', linetype = 'dashed')+ geom_rect(xmin = as.Date('2021-09-24'), xmax = as.Date('2021-12-17'), ymin=-Inf, ymax=Inf, alpha=0.005, fill="firebrick") + geom_rect(xmin = as.Date('2022-08-19'), xmax = as.Date('2022-11-25'), ymin=-Inf, ymax=Inf, alpha=0.005, fill="firebrick")+ scale_x_date(date_breaks = "1 month", date_labels = "%Y-%m")

ggsave(foreign_econ_dynamic_plot, file = here('figure_C4.pdf'), width = 7, height = 6)

#### Robustness: foreign econ wrt whole text (baseline = 2021 December) ----

coverage_foreignecon_dynamic <- feols(foreign_econ_news_share_wholetext ~ progov*i(month_yearDate, as.Date('2021-12-01')) | article_date, data=daily_trends_combined, panel.id = c('article_source', 'article_date'), vcov = 'DK')

coverage_foreignecon_dynamic_result <- summary(coverage_foreignecon_dynamic)

coverage_foreignecon_dynamic_result_data <- data.frame(coefs = coverage_foreignecon_dynamic_result$coefficients, se = coverage_foreignecon_dynamic_result$se) %>% tibble::rownames_to_column() %>% filter(grepl('progov:month_yearDate', rowname))


coverage_foreignecon_dynamic_result_data$lo95 <- coverage_foreignecon_dynamic_result_data$coefs - 1.96*coverage_foreignecon_dynamic_result_data$se 

coverage_foreignecon_dynamic_result_data$hi95 <- coverage_foreignecon_dynamic_result_data$coefs + 1.96*coverage_foreignecon_dynamic_result_data$se 

coverage_foreignecon_dynamic_result_data$variable_name <- stringr::str_split_fixed(coverage_foreignecon_dynamic_result_data$rowname, '::', 3)[, 2]

coverage_foreignecon_dynamic_result_data$variable_name <- as.Date(coverage_foreignecon_dynamic_result_data$variable_name)

coverage_foreignecon_dynamic_result_data <- coverage_foreignecon_dynamic_result_data %>% add_row(rowname = 'baseline', coefs = 0, se = 0, lo95 = 0, hi95 = 0, variable_name = as.Date('2021-12-01'))


foreign_econ_dynamic_plot <- ggplot(coverage_foreignecon_dynamic_result_data, aes(x=variable_name, y = coefs, ymin = lo95, ymax = hi95))+geom_point()+geom_linerange()+theme_bw()+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+ylab('Estimated effect of publishing an article on foreign economy')+ xlab('')+geom_hline(yintercept = 0, color = 'red', linetype = 'dashed')+ geom_rect(xmin = as.Date('2021-09-24'), xmax = as.Date('2021-12-17'), ymin=-Inf, ymax=Inf, alpha=0.005, fill="firebrick") + geom_rect(xmin = as.Date('2022-08-19'), xmax = as.Date('2022-11-25'), ymin=-Inf, ymax=Inf, alpha=0.005, fill="firebrick")+ scale_x_date(date_breaks = "1 month", date_labels = "%Y-%m")

ggsave(foreign_econ_dynamic_plot, file = here('figure_D1.pdf'), width = 7, height = 6)


# Exposed negativity analysis ----


#### Main: (baseline = 2021 December) ----

sentiment_dynamic <- feols(total_negative_exposure~ progov*i(month_yearDate, ref = as.Date('2021-12-01')) | article_date, data=daily_trends_combined, panel.id = c('article_source', 'article_date'), vcov = 'DK')

sentiment_dynamic_data <- data.frame(coefs = sentiment_dynamic$coefficients, se = sentiment_dynamic$se) %>% tibble::rownames_to_column() %>% filter(grepl('progov:month_yearDate', rowname))


sentiment_dynamic_data$lo95 <- sentiment_dynamic_data$coefs - 1.96*sentiment_dynamic_data$se 
sentiment_dynamic_data$hi95 <- sentiment_dynamic_data$coefs + 1.96*sentiment_dynamic_data$se 

sentiment_dynamic_data$variable_name <- stringr::str_split_fixed(sentiment_dynamic_data$rowname, '::', 3)[, 2]
sentiment_dynamic_data$variable_name <- as.Date(sentiment_dynamic_data$variable_name)

sentiment_dynamic_data <- sentiment_dynamic_data %>% add_row(rowname = 'baseline', coefs = 0, se = 0, lo95 = 0, hi95 = 0, variable_name = as.Date('2021-12-01'))


sentiment_dynamic_plot <- ggplot(sentiment_dynamic_data, aes(x=variable_name, y = coefs, ymin = lo95, ymax = hi95))+geom_point()+geom_linerange()+theme_bw()+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+ylab('Estimated exposed negativity in foreign economy news')+ xlab('')+geom_hline(yintercept = 0, color = 'red', linetype = 'dashed')+ geom_rect(xmin = as.Date('2021-09-24'), xmax = as.Date('2021-12-17'), ymin=-Inf, ymax=Inf, alpha=0.005, fill="firebrick") + geom_rect(xmin = as.Date('2022-08-19'), xmax = as.Date('2022-11-25'), ymin=-Inf, ymax=Inf, alpha=0.005, fill="firebrick")+ scale_x_date(date_breaks = "1 month", date_labels = "%Y-%m")


ggsave(sentiment_dynamic_plot, file = here('figure_5.pdf'), width = 7, height = 6)

# substantive interpretation
sentiment_dynamic_data[12:23, ] %>% pull(coefs) %>% mean()

#### Robustness: negativity wrt to the whole text(baseline = 2021 December) ----

sentiment_dynamic <- feols(total_negative_exposure_wholetext ~ progov*i(month_yearDate, ref = as.Date('2021-12-01')) | article_date, data=daily_trends_combined, panel.id = c('article_source', 'article_date'), vcov = 'DK')

sentiment_dynamic_data <- data.frame(coefs = sentiment_dynamic$coefficients, se = sentiment_dynamic$se) %>% tibble::rownames_to_column() %>% filter(grepl('progov:month_yearDate', rowname))


sentiment_dynamic_data$lo95 <- sentiment_dynamic_data$coefs - 1.96*sentiment_dynamic_data$se 
sentiment_dynamic_data$hi95 <- sentiment_dynamic_data$coefs + 1.96*sentiment_dynamic_data$se 

sentiment_dynamic_data$variable_name <- stringr::str_split_fixed(sentiment_dynamic_data$rowname, '::', 3)[, 2]
sentiment_dynamic_data$variable_name <- as.Date(sentiment_dynamic_data$variable_name)

sentiment_dynamic_data <- sentiment_dynamic_data %>% add_row(rowname = 'baseline', coefs = 0, se = 0, lo95 = 0, hi95 = 0, variable_name = as.Date('2021-12-01'))


sentiment_dynamic_plot <- ggplot(sentiment_dynamic_data, aes(x=variable_name, y = coefs, ymin = lo95, ymax = hi95))+geom_point()+geom_linerange()+theme_bw()+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+ylab('Estimated exposed negativity in foreign economy news')+ xlab('')+geom_hline(yintercept = 0, color = 'red', linetype = 'dashed')+ geom_rect(xmin = as.Date('2021-09-24'), xmax = as.Date('2021-12-17'), ymin=-Inf, ymax=Inf, alpha=0.005, fill="firebrick") + geom_rect(xmin = as.Date('2022-08-19'), xmax = as.Date('2022-11-25'), ymin=-Inf, ymax=Inf, alpha=0.005, fill="firebrick")+ scale_x_date(date_breaks = "1 month", date_labels = "%Y-%m")

ggsave(sentiment_dynamic_plot, file = here('figure_D3.pdf'), width = 7, height = 6)


#### Robustness: Exposed negativity (G7 only) ----

load(here('econnews_g7.rda'))


daily_sentiment_trends_g7 <- econ_news_g7 %>% group_by(article_source, article_date) %>% summarise(total_econ_news_negativity_g7 = sum(negative_econ_news, na.rm = TRUE))

daily_trends_combined <- full_join(daily_trends_combined, daily_sentiment_trends_g7, by = c('article_date', 'article_source'))

daily_trends_combined$total_econ_news_negativity_g7 <- ifelse(is.na(daily_trends_combined$total_econ_news_negativity_g7)==1, 0, daily_trends_combined$total_econ_news_negativity_g7)


daily_trends_combined$total_negative_exposure_g7 <- daily_trends_combined$total_econ_news_negativity_g7 / daily_trends_combined$total_daily_articles


sentiment_dynamic_month_g7 <- feols(total_negative_exposure_g7~ progov*i(month_yearDate, ref = as.Date('2021-12-01')) | article_date, data=daily_trends_combined, panel.id = c('article_source', 'article_date'), vcov = 'NW')


sentiment_dynamic_month_g7_result <- summary(sentiment_dynamic_month_g7)

sentiment_dynamic_data_month_g7 <- data.frame(coefs = sentiment_dynamic_month_g7_result$coefficients, se = sentiment_dynamic_month_g7_result$se) %>% tibble::rownames_to_column() %>% filter(grepl('progov:month_yearDate', rowname))

sentiment_dynamic_data_month_g7$lo95 <- sentiment_dynamic_data_month_g7$coefs - 1.96*sentiment_dynamic_data_month_g7$se 
sentiment_dynamic_data_month_g7$hi95 <- sentiment_dynamic_data_month_g7$coefs + 1.96*sentiment_dynamic_data_month_g7$se 

sentiment_dynamic_data_month_g7$variable_name <- stringr::str_split_fixed(sentiment_dynamic_data_month_g7$rowname, '::', 3)[, 2]
sentiment_dynamic_data_month_g7$variable_name <- as.Date(sentiment_dynamic_data_month_g7$variable_name)

sentiment_dynamic_data_month_g7 <- sentiment_dynamic_data_month_g7 %>% add_row(rowname = 'baseline', coefs = 0, se = 0, lo95 = 0, hi95 = 0, variable_name = as.Date('2021-12-01'))


econ_negativity_g7_plot <- ggplot(sentiment_dynamic_data_month_g7, aes(x=variable_name, y = coefs, ymin = lo95, ymax = hi95))+geom_point()+geom_linerange()+theme_bw()+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+ylab('Estimated exposed negativity in foreign economy news')+ xlab('')+geom_hline(yintercept = 0, color = 'red', linetype = 'dashed')+ geom_rect(xmin = as.Date('2021-09-24'), xmax = as.Date('2021-12-17'), ymin=-Inf, ymax=Inf, alpha=0.005, fill="firebrick") + geom_rect(xmin = as.Date('2022-08-19'), xmax = as.Date('2022-11-25'), ymin=-Inf, ymax=Inf, alpha=0.005, fill="firebrick")+ scale_x_date(date_breaks = "1 month", date_labels = "%Y-%m")

ggsave(econ_negativity_g7_plot, file = here('figure_E1.pdf'), width = 7, height = 6)


#### Robustness: Exposed negativity (Russia/Ukraine only) ----

load(here('econnews_russia_ukraine.rda'))


daily_sentiment_trends_russia_ukraine <- econ_news_Russia %>% group_by(article_source, article_date) %>% summarise(total_econ_news_negativity_ru = sum(negative_econ_news, na.rm = TRUE))

daily_trends_combined <- full_join(daily_trends_combined, daily_sentiment_trends_russia_ukraine, by = c('article_date', 'article_source'))

daily_trends_combined$total_econ_news_negativity_ru <- ifelse(is.na(daily_trends_combined$total_econ_news_negativity_ru)==1, 0, daily_trends_combined$total_econ_news_negativity_ru)


daily_trends_combined$total_negative_exposure_ru <- daily_trends_combined$total_econ_news_negativity_ru / daily_trends_combined$total_daily_articles


sentiment_dynamic_month_ru <- feols(total_negative_exposure_ru~ progov*i(month_yearDate, ref = as.Date('2021-12-01')) | article_date, data=daily_trends_combined, panel.id = c('article_source', 'article_date'), vcov = 'NW')


sentiment_dynamic_month_ru_result <- summary(sentiment_dynamic_month_ru)

sentiment_dynamic_data_month_ru <- data.frame(coefs = sentiment_dynamic_month_ru_result$coefficients, se = sentiment_dynamic_month_ru_result$se) %>% tibble::rownames_to_column() %>% filter(grepl('progov:month_yearDate', rowname))

sentiment_dynamic_data_month_ru$lo95 <- sentiment_dynamic_data_month_ru$coefs - 1.96*sentiment_dynamic_data_month_ru$se 
sentiment_dynamic_data_month_ru$hi95 <- sentiment_dynamic_data_month_ru$coefs + 1.96*sentiment_dynamic_data_month_ru$se 

sentiment_dynamic_data_month_ru$variable_name <- stringr::str_split_fixed(sentiment_dynamic_data_month_ru$rowname, '::', 3)[, 2]
sentiment_dynamic_data_month_ru$variable_name <- as.Date(sentiment_dynamic_data_month_ru$variable_name)

sentiment_dynamic_data_month_ru <- sentiment_dynamic_data_month_ru %>% add_row(rowname = 'baseline', coefs = 0, se = 0, lo95 = 0, hi95 = 0, variable_name = as.Date('2021-12-01'))


econ_negativity_ru_plot <- ggplot(sentiment_dynamic_data_month_ru, aes(x=variable_name, y = coefs, ymin = lo95, ymax = hi95))+geom_point()+geom_linerange()+theme_bw()+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+ylab('Estimated exposed negativity in foreign economy news')+ xlab('')+geom_hline(yintercept = 0, color = 'red', linetype = 'dashed')+ geom_rect(xmin = as.Date('2021-09-24'), xmax = as.Date('2021-12-17'), ymin=-Inf, ymax=Inf, alpha=0.005, fill="firebrick") + geom_rect(xmin = as.Date('2022-08-19'), xmax = as.Date('2022-11-25'), ymin=-Inf, ymax=Inf, alpha=0.005, fill="firebrick")+ scale_x_date(date_breaks = "1 month", date_labels = "%Y-%m")

ggsave(econ_negativity_ru_plot, file = here('figure_E2.pdf'), width = 7, height = 6)

#### Robustness: Exposed negativity (not G7 or Russia/Ukraine) ----

load(here('econnews_intl_others.rda'))


daily_sentiment_trends_intl_others <- econ_news_intl_others %>% group_by(article_source, article_date) %>% summarise(total_econ_news_negativity_intl_others = sum(negative_econ_news, na.rm = TRUE))

daily_trends_combined <- full_join(daily_trends_combined, daily_sentiment_trends_intl_others, by = c('article_date', 'article_source'))

daily_trends_combined$total_econ_news_negativity_intl_others <- ifelse(is.na(daily_trends_combined$total_econ_news_negativity_intl_others)==1, 0, daily_trends_combined$total_econ_news_negativity_intl_others)


daily_trends_combined$total_negative_exposure_intl_others <- daily_trends_combined$total_econ_news_negativity_intl_others / daily_trends_combined$total_daily_articles


sentiment_dynamic_month_intlothers <- feols(total_negative_exposure_intl_others~ progov*i(month_yearDate, ref = as.Date('2021-12-01')) | article_date, data=daily_trends_combined, panel.id = c('article_source', 'article_date'), vcov = 'NW')


sentiment_dynamic_month_intlothers_result <- summary(sentiment_dynamic_month_intlothers)

sentiment_dynamic_data_month_intlothers <- data.frame(coefs = sentiment_dynamic_month_intlothers_result$coefficients, se = sentiment_dynamic_month_intlothers_result$se) %>% tibble::rownames_to_column() %>% filter(grepl('progov:month_yearDate', rowname))

sentiment_dynamic_data_month_intlothers$lo95 <- sentiment_dynamic_data_month_intlothers$coefs - 1.96*sentiment_dynamic_data_month_intlothers$se 
sentiment_dynamic_data_month_intlothers$hi95 <- sentiment_dynamic_data_month_intlothers$coefs + 1.96*sentiment_dynamic_data_month_intlothers$se 

sentiment_dynamic_data_month_intlothers$variable_name <- stringr::str_split_fixed(sentiment_dynamic_data_month_intlothers$rowname, '::', 3)[, 2]
sentiment_dynamic_data_month_intlothers$variable_name <- as.Date(sentiment_dynamic_data_month_intlothers$variable_name)

sentiment_dynamic_data_month_intlothers <- sentiment_dynamic_data_month_intlothers %>% add_row(rowname = 'baseline', coefs = 0, se = 0, lo95 = 0, hi95 = 0, variable_name = as.Date('2021-12-01'))


econ_negativity_intlothers_plot <- ggplot(sentiment_dynamic_data_month_intlothers, aes(x=variable_name, y = coefs, ymin = lo95, ymax = hi95))+geom_point()+geom_linerange()+theme_bw()+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+ylab('Estimated exposed negativity in foreign economy news')+ xlab('')+geom_hline(yintercept = 0, color = 'red', linetype = 'dashed')+ geom_rect(xmin = as.Date('2021-09-24'), xmax = as.Date('2021-12-17'), ymin=-Inf, ymax=Inf, alpha=0.005, fill="firebrick") + geom_rect(xmin = as.Date('2022-08-19'), xmax = as.Date('2022-11-25'), ymin=-Inf, ymax=Inf, alpha=0.005, fill="firebrick")+ scale_x_date(date_breaks = "1 month", date_labels = "%Y-%m")

ggsave(econ_negativity_intlothers_plot, file = here('figure_E3.pdf'), width = 7, height = 6)



#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@


# World news coverage analysis ----

#### Main: (baseline = 2021 December) ----

world_dynamic <- feols(world_news_share~ progov*i(month_yearDate, ref = as.Date('2021-12-01')) | article_date, data=daily_trends_combined, panel.id = c('article_source', 'article_date'), vcov = 'DK')


world_dynamic_data <- data.frame(coefs = world_dynamic$coefficients, se = world_dynamic$se) %>% tibble::rownames_to_column() %>% filter(grepl('progov:month_yearDate', rowname))

world_dynamic_data$lo95 <- world_dynamic_data$coefs - 1.96*world_dynamic_data$se 
world_dynamic_data$hi95 <- world_dynamic_data$coefs + 1.96*world_dynamic_data$se 

world_dynamic_data$variable_name <- stringr::str_split_fixed(world_dynamic_data$rowname, '::', 3)[, 2]
world_dynamic_data$variable_name <- as.Date(world_dynamic_data$variable_name)

world_dynamic_data <- world_dynamic_data %>% add_row(rowname = 'baseline', coefs = 0, se = 0, lo95 = 0, hi95 = 0, variable_name = as.Date('2021-12-01'))


world_dynamic_plot <- ggplot(world_dynamic_data, aes(x=variable_name, y = coefs, ymin = lo95, ymax = hi95))+geom_point()+geom_linerange()+theme_bw()+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+ylab('Estimated effect of publishing an article on world news \n (excluding economy)')+ xlab('')+geom_hline(yintercept = 0, color = 'red', linetype = 'dashed')+ geom_rect(xmin = as.Date('2021-09-24'), xmax = as.Date('2021-12-17'), ymin=-Inf, ymax=Inf, alpha=0.005, fill="firebrick") + geom_rect(xmin = as.Date('2022-08-19'), xmax = as.Date('2022-11-25'), ymin=-Inf, ymax=Inf, alpha=0.005, fill="firebrick")+ scale_x_date(date_breaks = "1 month", date_labels = "%Y-%m")

ggsave(world_dynamic_plot, file = here('figure_B2.pdf'), width = 7, height = 6)


#### Robustness: (baseline = 2021 December) ----

world_dynamic <- feols(world_news_share_wholetext ~ progov*i(month_yearDate, ref = as.Date('2021-12-01')) | article_date, data=daily_trends_combined, panel.id = c('article_source', 'article_date'), vcov = 'DK')


world_dynamic_data <- data.frame(coefs = world_dynamic$coefficients, se = world_dynamic$se) %>% tibble::rownames_to_column() %>% filter(grepl('progov:month_yearDate', rowname))

world_dynamic_data$lo95 <- world_dynamic_data$coefs - 1.96*world_dynamic_data$se 
world_dynamic_data$hi95 <- world_dynamic_data$coefs + 1.96*world_dynamic_data$se 

world_dynamic_data$variable_name <- stringr::str_split_fixed(world_dynamic_data$rowname, '::', 3)[, 2]
world_dynamic_data$variable_name <- as.Date(world_dynamic_data$variable_name)

world_dynamic_data <- world_dynamic_data %>% add_row(rowname = 'baseline', coefs = 0, se = 0, lo95 = 0, hi95 = 0, variable_name = as.Date('2021-12-01'))


world_dynamic_plot <- ggplot(world_dynamic_data, aes(x=variable_name, y = coefs, ymin = lo95, ymax = hi95))+geom_point()+geom_linerange()+theme_bw()+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+ylab('Estimated effect of publishing an article on world news \n (excluding economy)')+ xlab('')+geom_hline(yintercept = 0, color = 'red', linetype = 'dashed')+ geom_rect(xmin = as.Date('2021-09-24'), xmax = as.Date('2021-12-17'), ymin=-Inf, ymax=Inf, alpha=0.005, fill="firebrick") + geom_rect(xmin = as.Date('2022-08-19'), xmax = as.Date('2022-11-25'), ymin=-Inf, ymax=Inf, alpha=0.005, fill="firebrick")+ scale_x_date(date_breaks = "1 month", date_labels = "%Y-%m")

ggsave(world_dynamic_plot, file = here('figure_D2.pdf'), width = 7, height = 6)


# Appendix: President Speech Analysis ----

load(here('pres_speech.rda'))

# Define time windows
date_windows <- list(
  c("2020-12-01", "2021-04-01"),
  c("2021-03-01", "2021-07-01"),
  c("2021-06-01", "2021-10-01"),
  c("2021-09-01", "2022-01-01"),
  c("2021-12-01", "2022-04-01"),
  c("2022-03-01", "2022-07-01"),
  c("2022-06-01", "2022-10-01"),
  c("2022-09-01", "2023-01-01"),
  c("2022-12-01", "2023-04-01")
)

# Filter and model for each window
models <- map(date_windows, ~ {
  df <- president_speech %>% filter(month_yearDate > .x[1], month_yearDate < .x[2], econ_statement == 1)
  
  feols(negative_econ_statement ~ foreign_mention + turkey_mention + domestic_company + domestic_facility, data = df)
})

# Compute average outcome for each window
avg_outcomes <- map_chr(date_windows, ~ {
  df <- president_speech %>%
    filter(
      month_yearDate > .x[1],
      month_yearDate < .x[2],
      econ_statement == 1
    )
  round(mean(df$negative_econ_statement, na.rm = TRUE), 2) %>% as.character()
})

# Create summary table
rows <- tibble(
  term = "Average outcome",
  !!!set_names(avg_outcomes, as.character(1:9))
)


names(models) <- c(
  "2021-Q1","2021-Q2","2021-Q3","2021-Q4",
  "2022-Q1","2022-Q2","2022-Q3","2022-Q4","2023-Q1"
)

modelsummary(models, stars_note = TRUE, stars = c('*' = .1, '**' = .05, '***' = 0.01), metrics = 'R2', gof_omit = 'R2 Within|R2 Adj|Std.Errors|FE', add_rows = rows, coef_map = c('foreign_mention' = 'Foreign nation/company mentioned', 'turkey_mention' = 'Turkey mentioned', 'domestic_company' = 'Domestic company mentioned', 'domestic_facility' = 'Domestic facility mentioned'), output = here('table_G1.txt'))


# Appendix: Chunk Analysis ----


### Domestic econ coverage after foreign economy coverage ----

domestic_coverage_m1 <- feols(domestic_follows_foreignecon_share~ progov*i(month_yearDate, ref = as.Date('2021-12-01')) | article_date, data=daily_trends_combined, panel.id = c('article_source', 'article_date'), vcov = 'DK')

domestic_coverage_m1_result <- summary(domestic_coverage_m1)

domestic_coverage_m1_result_data <- data.frame(coefs = domestic_coverage_m1_result$coefficients, se = domestic_coverage_m1_result$se) %>% tibble::rownames_to_column() %>% filter(grepl('progov:month_yearDate', rowname))


domestic_coverage_m1_result_data$lo95 <- domestic_coverage_m1_result_data$coefs - 1.96*domestic_coverage_m1_result_data$se 
domestic_coverage_m1_result_data$hi95 <- domestic_coverage_m1_result_data$coefs + 1.96*domestic_coverage_m1_result_data$se 

domestic_coverage_m1_result_data$variable_name <- stringr::str_split_fixed(domestic_coverage_m1_result_data$rowname, '::', 3)[, 2]
domestic_coverage_m1_result_data$variable_name <- as.Date(domestic_coverage_m1_result_data$variable_name)

domestic_coverage_m1_result_data <- domestic_coverage_m1_result_data %>% add_row(rowname = 'baseline', coefs = 0, se = 0, lo95 = 0, hi95 = 0, variable_name = as.Date('2021-12-01'))


domestic_coverage_plot <- ggplot(domestic_coverage_m1_result_data, aes(x=variable_name, y = coefs, ymin = lo95, ymax = hi95))+geom_point()+geom_linerange()+theme_bw()+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+ylab('Estimated effect of covering the domestic economy after \n covering a foreign economy \n in the same article')+ xlab('')+geom_hline(yintercept = 0, color = 'red', linetype = 'dashed')+ geom_rect(xmin = as.Date('2021-09-24'), xmax = as.Date('2021-12-17'), ymin=-Inf, ymax=Inf, alpha=0.005, fill="firebrick") + geom_rect(xmin = as.Date('2022-08-19'), xmax = as.Date('2022-11-25'), ymin=-Inf, ymax=Inf, alpha=0.005, fill="firebrick")+ scale_x_date(date_breaks = "1 month", date_labels = "%Y-%m")

## substantive interpretation:
round(mean(daily_trends_combined$domestic_follows_foreignecon_share), digits = 3)

round(domestic_coverage_m1_result_data[12:23, ] %>% pull(coefs) %>% mean(), digits = 3)


ggsave(domestic_coverage_plot, file = here('figure_H1.pdf'), width = 7, height = 6)



### Domestic neutral/positive econ coverage after negative foreign economy coverage ----

domestic_coverage_m2 <- feols(domestic_pos_neut_follows_foreign_neg_share~ progov*i(month_yearDate, ref = as.Date('2021-12-01')) | article_date, data=daily_trends_combined, panel.id = c('article_source', 'article_date'), vcov = 'DK')

domestic_coverage_m2_result <- summary(domestic_coverage_m2)

domestic_coverage_m2_result_data <- data.frame(coefs = domestic_coverage_m2_result$coefficients, se = domestic_coverage_m2_result$se) %>% tibble::rownames_to_column() %>% filter(grepl('progov:month_yearDate', rowname))


domestic_coverage_m2_result_data$lo95 <- domestic_coverage_m2_result_data$coefs - 1.96*domestic_coverage_m2_result_data$se 
domestic_coverage_m2_result_data$hi95 <- domestic_coverage_m2_result_data$coefs + 1.96*domestic_coverage_m2_result_data$se 

domestic_coverage_m2_result_data$variable_name <- stringr::str_split_fixed(domestic_coverage_m2_result_data$rowname, '::', 3)[, 2]
domestic_coverage_m2_result_data$variable_name <- as.Date(domestic_coverage_m2_result_data$variable_name)

domestic_coverage_m2_result_data <- domestic_coverage_m2_result_data %>% add_row(rowname = 'baseline', coefs = 0, se = 0, lo95 = 0, hi95 = 0, variable_name = as.Date('2021-12-01'))


domestic_neg_pos_coverage_plot <- ggplot(domestic_coverage_m2_result_data, aes(x=variable_name, y = coefs, ymin = lo95, ymax = hi95))+geom_point()+geom_linerange()+theme_bw()+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+ylab('Estimated effect of a neutral or positive domestic economy coverage \n after a negative foreign economy coverage \n in the same article')+ xlab('')+geom_hline(yintercept = 0, color = 'red', linetype = 'dashed')+ geom_rect(xmin = as.Date('2021-09-24'), xmax = as.Date('2021-12-17'), ymin=-Inf, ymax=Inf, alpha=0.005, fill="firebrick") + geom_rect(xmin = as.Date('2022-08-19'), xmax = as.Date('2022-11-25'), ymin=-Inf, ymax=Inf, alpha=0.005, fill="firebrick")+ scale_x_date(date_breaks = "1 month", date_labels = "%Y-%m")


ggsave(domestic_neg_pos_coverage_plot, file = here('figure_H2.pdf'), width = 7, height = 6)

#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

# Appendix: World news and foreign news coverage analysis ----

## G7 countries (econ news) ----

daily_g7_trends_econ <- econ_news_g7 %>% group_by(article_source, article_date) %>% summarise(total_daily_articles_g7_econ = n())

daily_trends_combined <- left_join(daily_trends_combined, daily_g7_trends_econ, by = c('article_source', 'article_date'))

daily_trends_combined$g7_share_econ <- daily_trends_combined$total_daily_articles_g7_econ / daily_trends_combined$total_daily_articles

daily_trends_combined$g7_share_econ <- ifelse(is.na(daily_trends_combined$g7_share_econ)==1, 0, daily_trends_combined$g7_share_econ)


coverage_g7 <- feols(g7_share_econ~ progov*i(month_yearDate, ref = as.Date('2021-12-01')) | article_date, data=daily_trends_combined, panel.id = c('article_source', 'article_date'), vcov = 'DK')

coverage_g7_result <- summary(coverage_g7)

coverage_data_g7 <- data.frame(coefs = coverage_g7_result$coefficients, se = coverage_g7_result$se) %>% tibble::rownames_to_column() %>% filter(grepl('progov:month_yearDate', rowname))

coverage_data_g7$lo95 <- coverage_data_g7$coefs - 1.96*coverage_data_g7$se 
coverage_data_g7$hi95 <- coverage_data_g7$coefs + 1.96*coverage_data_g7$se 

coverage_data_g7$variable_name <- stringr::str_split_fixed(coverage_data_g7$rowname, '::', 3)[, 2]
coverage_data_g7$variable_name <- as.Date(coverage_data_g7$variable_name)

coverage_data_g7 <- coverage_data_g7 %>% add_row(rowname = 'baseline', coefs = 0, se = 0, lo95 = 0, hi95 = 0, variable_name = as.Date('2021-12-01'))


round(coverage_data_g7[12:23, ] %>% pull(coefs) %>% mean(), digits = 2)
mean(daily_trends_combined$g7_share_econ)


coverage_plot_g7 <- ggplot(coverage_data_g7, aes(x=variable_name, y = coefs, ymin = lo95, ymax = hi95))+geom_point()+geom_linerange()+theme_bw()+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+ylab('Estimated effect of publishing economic news on G7 countries')+ xlab('')+geom_hline(yintercept = 0, color = 'red', linetype = 'dashed')+ geom_rect(xmin = as.Date('2021-09-24'), xmax = as.Date('2021-12-17'), ymin=-Inf, ymax=Inf, alpha=0.005, fill="firebrick") + geom_rect(xmin = as.Date('2022-08-19'), xmax = as.Date('2022-11-25'), ymin=-Inf, ymax=Inf, alpha=0.005, fill="firebrick")+ scale_x_date(date_breaks = "1 month", date_labels = "%Y-%m")

ggsave(coverage_plot_g7, file = here('figure_I1.pdf'), width = 7, height = 6)


## G7 countries (world news) ----

load(here('worldnews_g7.rda'))

# calculate counts each day-source and then merge with daily counts
daily_g7_trends <- world_news_g7 %>% group_by(article_source, article_date) %>% summarise(total_daily_articles_g7 = n())

daily_trends_combined <- left_join(daily_trends_combined, daily_g7_trends, by = c('article_source', 'article_date'))


daily_trends_combined$g7_share <- daily_trends_combined$total_daily_articles_g7 / daily_trends_combined$total_daily_articles

daily_trends_combined$g7_share <- ifelse(is.na(daily_trends_combined$g7_share)==1, 0, daily_trends_combined$g7_share)



coverage_g7 <- feols(g7_share~ progov*i(month_yearDate, ref = as.Date('2021-12-01')) | article_date, data=daily_trends_combined, panel.id = c('article_source', 'article_date'), vcov = 'DK')

coverage_g7_result <- summary(coverage_g7)

coverage_data_g7 <- data.frame(coefs = coverage_g7_result$coefficients, se = coverage_g7_result$se) %>% tibble::rownames_to_column() %>% filter(grepl('progov:month_yearDate', rowname))

coverage_data_g7$lo95 <- coverage_data_g7$coefs - 1.96*coverage_data_g7$se 
coverage_data_g7$hi95 <- coverage_data_g7$coefs + 1.96*coverage_data_g7$se 

coverage_data_g7$variable_name <- stringr::str_split_fixed(coverage_data_g7$rowname, '::', 3)[, 2]
coverage_data_g7$variable_name <- as.Date(coverage_data_g7$variable_name)

coverage_data_g7 <- coverage_data_g7 %>% add_row(rowname = 'baseline', coefs = 0, se = 0, lo95 = 0, hi95 = 0, variable_name = as.Date('2021-12-01'))


coverage_plot_g7 <- ggplot(coverage_data_g7, aes(x=variable_name, y = coefs, ymin = lo95, ymax = hi95))+geom_point()+geom_linerange()+theme_bw()+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+ylab('Estimated effect of publishing world news on G7 countries')+ xlab('')+geom_hline(yintercept = 0, color = 'red', linetype = 'dashed')+ geom_rect(xmin = as.Date('2021-09-24'), xmax = as.Date('2021-12-17'), ymin=-Inf, ymax=Inf, alpha=0.005, fill="firebrick") + geom_rect(xmin = as.Date('2022-08-19'), xmax = as.Date('2022-11-25'), ymin=-Inf, ymax=Inf, alpha=0.005, fill="firebrick")+ scale_x_date(date_breaks = "1 month", date_labels = "%Y-%m")

ggsave(coverage_plot_g7, file = here('figure_I2.pdf'), width = 7, height = 6)



#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@


## Russia/Ukraine news (econ news) ----

daily_ru_trends_econ <- econ_news_Russia %>% group_by(article_source, article_date) %>% summarise(total_daily_articles_ru_econ = n())

daily_trends_combined <- left_join(daily_trends_combined, daily_ru_trends_econ, by = c('article_source', 'article_date'))

daily_trends_combined$ru_share_econ <- daily_trends_combined$total_daily_articles_ru_econ / daily_trends_combined$total_daily_articles

daily_trends_combined$ru_share_econ <- ifelse(is.na(daily_trends_combined$ru_share_econ)==1, 0, daily_trends_combined$ru_share_econ)


coverage_ru <- feols(ru_share_econ~ progov*i(month_yearDate, ref = as.Date('2021-12-01')) | article_date, data=daily_trends_combined, panel.id = c('article_source', 'article_date'), vcov = 'DK')

coverage_ru_result <- summary(coverage_ru)

coverage_data_ru <- data.frame(coefs = coverage_ru_result$coefficients, se = coverage_ru_result$se) %>% tibble::rownames_to_column() %>% filter(grepl('progov:month_yearDate', rowname))

coverage_data_ru$lo95 <- coverage_data_ru$coefs - 1.96*coverage_data_ru$se 
coverage_data_ru$hi95 <- coverage_data_ru$coefs + 1.96*coverage_data_ru$se 

coverage_data_ru$variable_name <- stringr::str_split_fixed(coverage_data_ru$rowname, '::', 3)[, 2]
coverage_data_ru$variable_name <- as.Date(coverage_data_ru$variable_name)

coverage_data_ru <- coverage_data_ru %>% add_row(rowname = 'baseline', coefs = 0, se = 0, lo95 = 0, hi95 = 0, variable_name = as.Date('2021-12-01'))


round(coverage_data_ru[12:23, ] %>% pull(coefs) %>% mean(), digits = 2)
mean(daily_trends_combined$ru_share_econ)


coverage_plot_ru <- ggplot(coverage_data_ru, aes(x=variable_name, y = coefs, ymin = lo95, ymax = hi95))+geom_point()+geom_linerange()+theme_bw()+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+ylab('Estimated effect of publishing economic news on Russia/Ukraine')+ xlab('')+geom_hline(yintercept = 0, color = 'red', linetype = 'dashed')+ geom_rect(xmin = as.Date('2021-09-24'), xmax = as.Date('2021-12-17'), ymin=-Inf, ymax=Inf, alpha=0.005, fill="firebrick") + geom_rect(xmin = as.Date('2022-08-19'), xmax = as.Date('2022-11-25'), ymin=-Inf, ymax=Inf, alpha=0.005, fill="firebrick")+ scale_x_date(date_breaks = "1 month", date_labels = "%Y-%m")

ggsave(coverage_plot_ru, file = here('figure_I3.pdf'), width = 7, height = 6)


## Russia/Ukraine news (world news) ----

load(here('worldnews_russia_ukraine.rda'))

# calculate counts each day-source and then merge with daily counts
daily_ru_trends <- world_news_Russia %>% group_by(article_source, article_date) %>% summarise(total_daily_articles_ru = n())

daily_trends_combined <- left_join(daily_trends_combined, daily_ru_trends, by = c('article_source', 'article_date'))


daily_trends_combined$ru_share <- daily_trends_combined$total_daily_articles_ru / daily_trends_combined$total_daily_articles

daily_trends_combined$ru_share <- ifelse(is.na(daily_trends_combined$ru_share)==1, 0, daily_trends_combined$ru_share)



coverage_ru <- feols(ru_share~ progov*i(month_yearDate, ref = as.Date('2021-12-01')) | article_date, data=daily_trends_combined, panel.id = c('article_source', 'article_date'), vcov = 'DK')

coverage_ru_result <- summary(coverage_ru)

coverage_data_ru <- data.frame(coefs = coverage_ru_result$coefficients, se = coverage_ru_result$se) %>% tibble::rownames_to_column() %>% filter(grepl('progov:month_yearDate', rowname))

coverage_data_ru$lo95 <- coverage_data_ru$coefs - 1.96*coverage_data_ru$se 
coverage_data_ru$hi95 <- coverage_data_ru$coefs + 1.96*coverage_data_ru$se 

coverage_data_ru$variable_name <- stringr::str_split_fixed(coverage_data_ru$rowname, '::', 3)[, 2]
coverage_data_ru$variable_name <- as.Date(coverage_data_ru$variable_name)

coverage_data_ru <- coverage_data_ru %>% add_row(rowname = 'baseline', coefs = 0, se = 0, lo95 = 0, hi95 = 0, variable_name = as.Date('2021-12-01'))


coverage_plot_ru <- ggplot(coverage_data_ru, aes(x=variable_name, y = coefs, ymin = lo95, ymax = hi95))+geom_point()+geom_linerange()+theme_bw()+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+ylab('Estimated effect of publishing world news on Russia/Ukraine')+ xlab('')+geom_hline(yintercept = 0, color = 'red', linetype = 'dashed')+ geom_rect(xmin = as.Date('2021-09-24'), xmax = as.Date('2021-12-17'), ymin=-Inf, ymax=Inf, alpha=0.005, fill="firebrick") + geom_rect(xmin = as.Date('2022-08-19'), xmax = as.Date('2022-11-25'), ymin=-Inf, ymax=Inf, alpha=0.005, fill="firebrick")+ scale_x_date(date_breaks = "1 month", date_labels = "%Y-%m")

ggsave(coverage_plot_ru, file = here('figure_I4.pdf'), width = 7, height = 6)


## War\conflict analysis ----

load(here('war_news.rda'))

daily_trends_war <- war_news %>% group_by(article_source, article_date) %>% summarise(total_daily_articles_war = n())

daily_trends_combined <- left_join(daily_trends_combined, daily_trends_war, by = c('article_source', 'article_date'))


daily_trends_combined$war_share <- daily_trends_combined$total_daily_articles_war / daily_trends_combined$total_daily_articles

daily_trends_combined$war_share <- ifelse(is.na(daily_trends_combined$war_share)==1, 0, daily_trends_combined$war_share)


war_model <- feols(war_share ~ progov*i(month_yearDate, ref = as.Date('2021-12-01')) | article_date, data=daily_trends_combined, panel.id = c('article_source', 'article_date'), vcov = 'DK')

war_model_result <- summary(war_model)

war_model_data <- data.frame(coefs = war_model_result$coefficients, se = war_model_result$se) %>% tibble::rownames_to_column() %>% filter(grepl('progov:month_yearDate', rowname))

war_model_data$lo95 <- war_model_data$coefs - 1.96*war_model_data$se 
war_model_data$hi95 <- war_model_data$coefs + 1.96*war_model_data$se 

war_model_data$variable_name <- stringr::str_split_fixed(war_model_data$rowname, '::', 3)[, 2]
war_model_data$variable_name <- as.Date(war_model_data$variable_name)

war_model_data <- war_model_data %>% add_row(rowname = 'baseline', coefs = 0, se = 0, lo95 = 0, hi95 = 0, variable_name = as.Date('2021-12-01'))


plot_war <- ggplot(war_model_data, aes(x=variable_name, y = coefs, ymin = lo95, ymax = hi95))+geom_point()+geom_linerange()+theme_bw()+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+ylab('Estimated effect of publishing news on war')+ xlab('')+geom_hline(yintercept = 0, color = 'red', linetype = 'dashed')+ geom_rect(xmin = as.Date('2021-09-24'), xmax = as.Date('2021-12-17'), ymin=-Inf, ymax=Inf, alpha=0.005, fill="firebrick") + geom_rect(xmin = as.Date('2022-08-19'), xmax = as.Date('2022-11-25'), ymin=-Inf, ymax=Inf, alpha=0.005, fill="firebrick")+ scale_x_date(date_breaks = "1 month", date_labels = "%Y-%m")

ggsave(plot_war, file = here('figure_I5.pdf'), width = 7, height = 6)



# Appendix: Survey Data ----


#### 2023 Survey Data ----
load(here('survey2023.rda'))

sabah_media_supporters <- d2023 %>% filter(sabah_media==1)

# the share of AKP-MHP voters among sabah media consumers:
sabah_progov <- mean(sabah_media_supporters$progov)


trt_media_supporters <- d2023 %>% filter(trt==1)

# the share of AKP-MHP voters among trt media consumers:
trt_progov <- mean(trt_media_supporters$progov)



sozcu_media_supporters <- d2023 %>% filter(sozcu_media==1)

# the share of AKP-MHP voters among sozcu consumers:
sozcu_progov <- mean(sozcu_media_supporters$progov)



party_media2023 <- data.frame(media = c('Sozcu', 'Sozcu', 'Sabah', 'Sabah', 'TRT', 'TRT'), party = c('Incumbent voter share', 'Non-incumbent voter share', 'Incumbent voter share', 'Non-incumbent voter share', 'Incumbent voter share', 'Non-incumbent voter share'), vote_share = c(sozcu_progov, 1-sozcu_progov, sabah_progov, 1-sabah_progov, trt_progov, 1-trt_progov))

party_media_plot <- ggplot(party_media2023, aes(x = media, y = vote_share, fill = party)) + geom_col(position = 'dodge') + theme_bw() + theme(legend.position = 'bottom') + labs(y = 'Composition of voters within each media outlet', x ='' , fill = '') + scale_fill_manual(values = c("Incumbent voter share" = "firebrick", "Non-incumbent voter share" = "dodgerblue4"))


ggsave(party_media_plot, file = here('figure_I7.pdf'), width = 6, height = 6)


issue_party_data1 <- d2023 %>% group_by(progov) %>% summarise(issue_prop = mean(econ_imp_combined))

issue_party_data1$issue <- 'Economy'

issue_party_data2 <- d2023 %>% group_by(progov) %>% summarise(issue_prop = mean(intl_imp_combined))

issue_party_data2$issue <- 'International affairs'

issue_party_data <- rbind(issue_party_data1, issue_party_data2)

issue_party_data$party <- ifelse(issue_party_data$progov == 1, 'Incumbent voters', 'Non-incumbent voters')

issue_party_plot <- ggplot(issue_party_data, aes(x = party, y = issue_prop, fill = issue)) + geom_col(position = 'dodge') + theme_bw() + theme(legend.position = 'bottom') + labs(y = 'Most/second most important issue in the country', x ='' , fill = '') + scale_fill_manual(values = c("Economy" = "firebrick", "International affairs" = "dodgerblue4"))

ggsave(issue_party_plot, file = here('figure_I9.pdf'), width = 6, height = 6)

#### 2018 Survey Data ----

load(here('survey2018.rda'))

sabah_media_supporters <- d2018 %>% filter(sabah_media==1)

# the share of AKP-MHP voters among sabah media consumers:
sabah_progov <- mean(sabah_media_supporters$progov)


trt_media_supporters <- d2018 %>% filter(trt==1)

# the share of AKP-MHP voters among trt media consumers:
trt_progov <- mean(trt_media_supporters$progov)



sozcu_media_supporters <- d2018 %>% filter(sozcu==1)

# the share of AKP-MHP voters among sozcu consumers:
sozcu_progov <- mean(sozcu_media_supporters$progov)



party_media2018 <- data.frame(media = c('Sozcu', 'Sozcu', 'Sabah', 'Sabah', 'TRT', 'TRT'), party = c('Incumbent voter share', 'Non-incumbent voter share', 'Incumbent voter share', 'Non-incumbent voter share', 'Incumbent voter share', 'Non-incumbent voter share'), vote_share = c(sozcu_progov, 1-sozcu_progov, sabah_progov, 1-sabah_progov, trt_progov, 1-trt_progov))

party_media_plot <- ggplot(party_media2018, aes(x = media, y = vote_share, fill = party)) + geom_col(position = 'dodge') + theme_bw() + theme(legend.position = 'bottom') + labs(y = 'Composition of voters within each media outlet', x ='' , fill = '') + scale_fill_manual(values = c("Incumbent voter share" = "firebrick", "Non-incumbent voter share" = "dodgerblue4"))


ggsave(party_media_plot, file = here('figure_I6.pdf'), width = 6, height = 6)

issue_party_data1 <- d2018 %>% group_by(progov) %>% summarise(issue_prop = mean(econ_imp_combined))

issue_party_data1$issue <- 'Economy'

issue_party_data2 <- d2018 %>% group_by(progov) %>% summarise(issue_prop = mean(intl_imp_combined))

issue_party_data2$issue <- 'International affairs'

issue_party_data <- rbind(issue_party_data1, issue_party_data2)

issue_party_data$party <- ifelse(issue_party_data$progov == 1, 'Incumbent voters', 'Non-incumbent voters')

issue_party_plot <- ggplot(issue_party_data, aes(x = party, y = issue_prop, fill = issue)) + geom_col(position = 'dodge') + theme_bw() + theme(legend.position = 'bottom') + labs(y = 'Most/second most important issue in the country', x ='' , fill = '') + scale_fill_manual(values = c("Economy" = "firebrick", "International affairs" = "dodgerblue4"))

ggsave(issue_party_plot, file = here('figure_I8.pdf'), width = 6, height = 6)

# Appendix: Generalized additive mixed model ----

daily_trends_combined <- daily_trends_combined %>% group_by(month_yearDate) %>%  mutate(month_id = cur_group_id()) 

# 0 won't be accepted in the function below, so re-create this:
daily_trends_combined$progov2 <- ifelse(daily_trends_combined$progov == 0, 1, 2)


daily_trends_combined$week_day <- wday(daily_trends_combined$article_date)

generalized_model <- gamm(foreign_econ_news_share~ s(month_id, by = progov2) + factor(week_day), data=daily_trends_combined, correlation = corAR1(form= ~article_date | article_source))

predictionData <- data.frame(month_id=c(rep(seq(1:26), 2)), progov2 = c(rep(1, 26), rep(2, 26)), week_day = 3)

genmodel_predictions <- predict(generalized_model$gam, newdata = predictionData, se.fit = TRUE, type = 'response') %>% data.frame()

genmodel_predictions <- cbind(genmodel_predictions, predictionData)

genmodel_predictions$lo95 <- genmodel_predictions$fit - 1.96*genmodel_predictions$se.fit

genmodel_predictions$hi95 <- genmodel_predictions$fit + 1.96*genmodel_predictions$se.fit

## Merge back month info
month_info <- daily_trends_combined %>% select(month_yearDate, month_id) %>% unique()

genmodel_predictions <- left_join(genmodel_predictions, month_info, by = 'month_id')


genmodel_predictions$progov_type <- ifelse(genmodel_predictions$progov2 == 1, 'Anti-gov outlet', 'Pro-gov outlet')

genmodel_plot <- ggplot(data = genmodel_predictions, aes(x = month_yearDate, y = fit, ymin = lo95, ymax = hi95, fill = factor(progov_type), color = factor(progov_type)))+geom_line()+geom_ribbon(alpha = 0.3)+theme_bw()+ geom_hline(yintercept = 0, color = 'black', linetype = 'dashed')+ scale_x_date(date_breaks = "1 month", date_labels = "%Y-%m") + theme(panel.grid.major.y = element_blank(), axis.text= element_text(size=18), axis.title=element_text(size=18), axis.text.x = element_text(angle = 45, hjust = 1, size = 14), legend.position = 'bottom')+ylab('Predicted probability of \n publishing an article on foreign economy')+xlab('') + labs(fill = '', color = '')


ggsave(genmodel_plot, file = here('figure_J1.pdf'), width = 8, height = 7)



